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Bayesian Analysis for Linearized 
Multi-Stage Models in Quantal Bioassay 

LYNN KUO and MICHAEL P. COHEN * 

Abstract 

Bayesian methods for estimating dose response curves in quantal bioas- 
say are studied. A linearized multi-stage model is assumed for the 
shape of the curves. A Gibbs sampling approach with data augmen- 
tation is employed to compute the Bayes estimates. In addition, esti- 
mation of the "relative additional risk" and the "risk specific dose" 
is studied. Model selection based on conditional predictive ordinates 
from cross- validated data is developed. 

KEY WORDS: Dose response curve; Gibbs sampling; Model selec- 
tion; Risk; Stochastic substitution. 
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!• Introduction 

Several models are available to explain the relationship between the re- 
sponse probability and the dosage in quantal bioassay. These include the 
probit, logit, one-hit, multi-hit, multistag^s and Weibull models. Krewski 
and Van Ryzin (1981) provide a review. The maximum likelihood method is 
usually employed to estimate the response curves. The Bayesian approach 
to these problems has been limited to a few simplified cases. Petrasovi- 
tas and Cornell (1975) propose approximations to the Bayes estimates for 
the one-hit model with single-dose assay. Zellner and Rossi (1982) employ 
numerical integration techniques for computing Bayes estimates in probit 
and logit models with noninformative priors. Recently, some advances have 
been made for more complex situations. Albert and Chib (1991) develop 
a Gibbs sampling approach for the probit model. They also extend their 
methods to polytomous responses and the i-link distribution. Gelfand and 
Kuo (1991) develop Gibbs sampling approaches for nonparametric dose 
response curves with Dirichlet process priors and with product of beta pri- 
ors. They also extend their models and methods to polytomous responses. 
Ramgopal, Laud, and Smith (1992) extend the work of Gelfand and Kuo 
to convex, concave, or ogive potency curves. Nagaraj (1991) proposes an 
importance sampling approach for the linearized multi-stage model. 

The linearized multi-stage model specifies the probability of a subject's 
positive response to a dose at d as 

P(cii/3) - 1 e^^+^»'+* '^^*'') (1.1) 

with Pj > 0, j = 0, This model is currently used by the U.S. 

Environmental Protection Agency (EPA) for risk assessment. As discussed 
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by Crump, Hoel, Langley, and Peto (1976) and Crump (1979), the model 
is based on the assumption that inductior of toxic effects results from a 
number of different random biological events, with the age-specific rate of 
occurrence of each event linearly related to dose. When A; = 1, this is the 
on^-hit model described by Armitage and Doll (1961). The model in (1.1) 
can also be considered as a generalized linear model (McCullagh and Nelder 
1989) with complementary log as the link fimction. 

To apply Bayesian analysis for the unknown parameters /3 = (/5o, . . . , ^/t)) 
we assume that the • - • 0k are independent, where fij has a ganuna dis- 
tribution r(aj,7j) with mean aj/jj for j = 0, . . . ,A;. The gamma distri- 
butions are chosen partly for convenience. The family of two parameter 
gamma densities, though, incorporates densities with both increasing and 
decreasing failure rates. Often statisticians have some prior belief about 
the location and spread of the distribution of 0j for each j. Let fij and cj 
denote a prior guess of this location and spread. Thei) aj and 7j can be 
chosen by the method of matching moments so that they satisfy fij = Oij/^j 
and (r^2 = q^/tJ. 

In risk assessment, subjects are divided into groups, say I groups, 
where for t = 1, . . . , /, rii subjects are given the substance at dosage level 
d,. The number of casualties observed in group i is denoted by x^. Let 
X = (xj, . . . ,x/). Given /3, the subjects' responses are assumed to be inde- 
pendent. Therefore the likeUhood function is 



It can be shown that the posterior distribution is a mixture of prod- 
ucts of gamma distributions. This mixture increases in complexity as / 
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increases. This paper provides a Gibbs sampling approach to computing 
the Bayes estimates. It can be used for any number of groups. This ap- 
proach augments the data by latent variables that convert the integrals 
(cumulative probabilities and survival probabilities) in the likelihood to 
probability density fimctions. This augmentation facilitates us in speci- 
fying the conditional density of (3 given the latent variables. A repeated 
sampling scheme that generates variates from this conditional density and 
the conditional density of the latent variables given /3 and the data allows 
us to approximate the posterior distribution. 

The dose response model in (1.1) is often used for low-dose extrapola- 
tion, where the potency of a toxic substance on human beings (at low dose) 
is inferred from experiments conducted at much higher doses. Practition- 
ers are often interested in "Relative Additional Risk" (RAR) and "Risk 
Specific Dose" (RSD). We provide methods of Bayesian inference for these 
quantities. RAR can be interpreted as the conditional probability of a sub- 
ject acquiring cancer, say, at dose level d, given that no cancer has been 
developed in the absence of the carcinogen. RSD measures the dosage such 
that the probability of a casualty is no more than a given number. 

Because the model in (1.1) depends on k, the number of parameters, 
it is natural to ask how to select k. One route to selecting k has been for 
biologists to propose theory to explain the model. Another route, more 
empiric al, is to fit the models to data. Chi-squared goodness of fit tests are 
usually employed by frequentists for model checking. A criterion for model 
selection based on Bayesian predictive ordinates from cross-validated data 
is adopted in this paper. 

Section 2 describes the Gibbs sampling approach. Section 3 addresses 
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estimation of RAR and RSD. Section 4 develops the criterion for model 
selection. Examples are given in Section 5, and some concluding remarks 
are provided in Section 6. 

2. Gibbs Sampling 

A Gibbs sampling approach specific for our problem is developed here 
to compute the posterior distribution of /3. One can refer to Geman and 
Geman (1984), Tanner and Wong (1987), and Gelfand and Smith (1990) 
for general discussion of Gibbs sampling. 

Given the likelihood and the prior, we observe that the posterior density 
of (3 is 

/(/3|x) oc n (l - e-^^'-'^'^^-^^^''^y n /?J^-^e-l^^+Sl.{-.-0^1^.. (2.1) 
t=l j=o 

To sample f3 from (2.1) by stochastic substitution with augmented data 
z, we can sample repeatedly from /(z|/3,x), the conditional distribution 
of z given (3 and x, and /(/3|z,x), the conditional distribution of /3 given 
z and X. Therefore, it is desirable to augment the data by the latent 
variables z, sucl) that we can easily sample variates from the two conditional 
densities /(z|Ax) and /(/3|z,x). For each i, let us consider the latent 
variables {zij}^ j = 1, . . . as i.i.d. random variables with an exponential 
distribution with mean l/(/3o + + Observe that 

P(z,,<l) = l-e-(^+^"'*+-+^*'f). 

A.ssume the distribution of the z,- = (z.i, . . • ,2.,i,) are independent for i = 
1, . . . , /. Let z = (zi, . . . ,z;). Then the conditional density of z given /3 
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and X is 

5(zlAx) = nn5i(^ol/3,x), (2.2) 

where 5, is derived from an exponential distribution with mean l/(/3o + 

fi^di + h/Sfcrff ) by truncation on the right by 1. To obtain the conditional 

density of (3 given z and x, set s,Q3) = /3o + M + ■ ■ ■ + Mi and observe 
/(/31z,x) « Y[/3^^-'e-bi+ZL(-<-'<m 



xn[s,(/3)]^'e-''-^&^"\ (2.3) 



.=1 



This suggests another data augmentation step is needed to facilitate the 
generation of the /3 variate given the latent variables. For each t, consider 
Wi = (wio, . . . ,Wik) distributed as a multinomial distribution with parame- 
ters X, and cell probabilities Moreover, generate 
the -wi, . . . , w; independently. Let w denote (wi, . . . , w;). Then the joint 
density of the w is given by 



; 

= n 



(2.4) 



Let z,+ = EjLi^ij- Now the conditional density of (3 given z,w, and x 
can be verified to be the product of independent gamma densities with 
parameters aj = ELi '^ij and ■yj = ELi di{zi+ + ni- x<), 

/031z,w,x)«n/5r'^~''^'- (2.5) 
Having explained the conditional densities needed in the substitution algo- 
ritlim, we list the steps used to generate the samplers. 
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(1) Start with an initial dioice of the A denoted by (^^\ Set the counter 
in the superscript to 0. 

(2) Generate z^'^ and w^'^ independently from /(z|/3^°\x) as in (2.2) and 
/(wl/3^°\x) as in (2.4). 

(3) Generate (^^^ from /(/3|z('), w(^),x) as in (2,5). 

(4) Repeat steps (1), (2), and (3) with the updated counter in all the 
superscripts, that is, with ^'^'^ replacing etc., in each step. 

The generation of z.j in (3) from the truncated density can be facilitated 
by use of the equation z.,- = G-'{G{0) + Uij{G{\)-G{0) )), where the Uij are 
independent uniform [0,1] variates, and G is the exponential distribution 
function with mean l/s,(/3). Then, z^,- = - ln(l - C/i;(l - e-"^))/si(/3), 
independently for j = i, . . . ,Xi. 

Suppose we also replicate the above iterations R times by using inde- 
pendent initial choices of the /3. Let the A z, and w variates generated in 
the L"" step of the iteration for the r"" replication be denoted by /3^/\z[^\ 
and for r = 1,...,/?.. Similar notation is applied to the variates 

Wtj,z,j.z,+ , etc. Then the posterior distribution of /3 is estimated by 

/03ix)« ^ Y:m4'^M'\^\ (2-6) 

^ r=l 

where / on the right is the product of gamma distributions described as 
in (2.5). Since the Pj, j = 0, . . . , A;, are independent in (2.5), the posterior 
mean of Pj can be approximated by 

.if ___^LtSk!fi . (2.7) 
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The computations in (2,6) and (2,7) use only the variates generated 
in the last step of each iteration. A more efficient algorithm has been 
proposed by Gelman and Rubin (1992) where estimates are produced by 
averaging over all variates generated in the second halves of the iterations 
for all replications. They also study the choices of L and R by analysis of 
variance techniques. 

We will use the empirical values {/3j-fr^}^=i from the repUcated samples 
to construct the predictive interval for 

3. Relative Additional Risk and Risk Specific 
Dose 

The computation of some functionals of the parameters by stochastic 
substitution c'lscussed in this section. 

The relative additional risk over the background at a dose of level d is 
defined by 

The relative additional risk can be interpreted as the conditional probability 
of acquiring cancer, for example, at dose level d given that no cancer has 
developed in the absence of the carcinogen. Observe that A{d) « at 
low dose. Let i^^^rf) denote an estimate of A{d) computed from the L"* 
iteration and the r*'' replication in the sampling procedure. Then A[^\d) 
can be computed from (2,5) and the moment generating functions of the 
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gaiiuna distributions 

/=0 V + 7; + E.il 4(.^i + ^.+!r") 



The Bayes estimate for >i(ci), denoted by i(d), is computed from 

Md)= ^ E^^^d)- (3-2) 

^ r=l 



The predictive interval for A{d) can be computed from the empirical mea- 
sure that assigns weight 1/R to eadi Ai^\d) for r = 1, . . . 

The Risk Specific Dose (RSD) is defined as the dose level such that the 
probability of a casualty P{d) is at most 9. RSD{9) is the -value d such 
that P{d) = 9. Equivalently, RSD{e) is defined to be the largest real root 
of the following equation: 

ln(l -0) + i3o + (3id+--- + M"" = 0. (3.3) 

Since the parameters {/3j}5=o are unknown, we propose to estimate the 
j = 0,...,/t, by (2.6) first, then find the roots of (3.3) where is 

replaced by /3j, the e.stimates of /3j, for all j = Q k. Observe that when 

k = l, then RSD{9) = -[/3o + ln(l - 

4. Model Selection 

The selection of A;, the degree of the polynomial in the exponent of (1.1), 
is studied in this section. A criterion based on the conditional predicted 
ordinates from cross-validated data is adopted for model selection. This 
approacii, called pseudo-marginal likelihood, was first proposed by Geisser 
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and Eddy (1979). The role of using predictive densities from cross-validated 
data for model determination was also supported by Stone (1974), Box 
(1980), and Pettit and Young (1990). Box (1980) states that the predictive 
distribution enables one to criticize the entertained model in the light of 
current data. Gelfand, Dey, and Chang (1992) have provided more detailed 
discussion of Bayesian model selection and its implementation by means of 
Gibbs sampling. 

Because the selection criterion was developed for independent observa- 
tions, we need to reconsider our data from this viewpoint. Suppose each 
subject is given a dose independently. The response of each subject is ei- 
ther positive or negative. Therefore, we observe n independent Bernoulli 
random variables with probability of success p^, s = 1, . . . , n. If d^s) is the 
dosage for subject 5, then 

Let y(5) denote the whole data set with the 5^'* observation deleted. The 
conditional predicted ordinate (CPO) for y,, denoted by c,, is defined by 

cs - /(y.iy(-)) / fiyMyis))f(P\yis))df3. (^-i) 

Note Cs depends on A:. We suppress this dependence in the notation for tlie 
time being for simphcity. 

The estimated density f(j3\y{s)) can be computed as in (2.6) with j/, 
deleted from the data set. Let q5J;J(,) and 7j>J(,) denote the parameters of 
the gamma density of (3j in (2.5) and (2.6) computed from the data with 
J/, deleted. Since /(j/.|Ay(5)) = /(y.l/3) = pK'(1 ~ P.)^^^% the computation 
of (4.1) is straightforward with the use of the moment generating function 
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of a gaimna distribution. If is 0 and the 6^^ subject is given the dosage 
d, then c, is approximated by 

,(L) 

If y, is 1 and the s"* subject is given the dosage d, then c, is approximated 
by one minus the right side of(4.2). The variance of CPO c, is given by 

Vic) = / e-2(^-^"'+ -^^"^)/(/3|y(.))d/3- (/ e-(^+^"'-^-^^'^*V(/3|y(.))ci/3) . 
It is approximated by 

.(/.) 

^- r=l j=0 \^ 7j,r,(3)' ' + ^(i^ j 

Since c, depends on k, we now use the notation c,{k) for r,, based on the 
model witli a k degree polynomial in the generalized linear model. Define 
the p.seudo-marginal likelihood C{k) = n','=i c,{k). We select the A; with the 
largest pseudo-marginal likelihood. Equivalently, we select the model with 
the largest Ey=i In c,{k). 

5. Numerical Examples 

5.1 Bayes Inference 

The data set is taken from Howe, Crump, and Van Landingham (1986). 
Nagaraj (1991) also used it for illustration. There are 200 animals divided 
into four groups with 50 animals in each group. The dosages assigned to the 
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groups are 0, 1, 2, and 3. The observed number of positive responses x in 
each group is (0,2,10,30). Tables 1-3 exhibit the various Bayes estimates 
computed based on the Gibbs sampUng method with 50 iterations and 1000 
replications. Three models with A: = 1,2,3 are considered, each with three 
choices of priors. The choices are (1) = 2 and 7j = 20 for all j, (2) 
a. =y. = 1 for all j, and (3) a_, = 1 and 7j = .01. The priors are chosen 
to center at .1, 1, and 100 with small, medium, and large variances. The 
parameters and functionals of interest are , j = 0, . . . , A;, and the Relative 
Additional Risk (RAR) Aid) evaluated at d = .5,1,2, and 3 (A(0) = 0 by 
definition). For eadi of these functionals, we display the point estimate 
(posterior mean), the posterior standard deviation, and the posterior 95% 
confidence interval. The Bayes estimates for P{d) evaluated at the dosage 
levels are also given. The corresponding standard deviations and confidence 
intervals can be computed from the Gibbs samplers by methods similar to 
those for /?_,■ and A{d). They are omitted here for brevity. The risk specific 
dose (RSD) for 9 = .025, .05, .075, .1, .125, and .15 is provided. 

The results of the.se tables show that the Bayes estimates for the second 
and third choices of the priors are very close to each other. The effect of 
the small mean and variance of the first prior can be seen from the highest 
coefEci^ at (/?*:) for each model. As k increases, the Bayes estimates of P(d) 
for d = 0,1,2, and 3 approach the nonparametric maximum likelihood 
estimates, which are 0, .04, .2, and .6, respectively. 

5.2 Model Selection 

The numerical results for the model selection are given here. We have 
considered each y,, s = 1, ...,200, to have a Bernoulli distribution with 
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the probability of a success (positive response) to be P(d), where d is the 
dose given to the subject s. Given the data x = (0,2,10,30) observed at 
the dosages 0, 1, 2, and 3, we can consider = 0 for 5 = 1-50, 53-100, 
111-150, 181-200, and y, = 1 for the rest of the indices. 

Because of the grouped data x, many of the CPO's will have the same 
value. In fact, instead of computing 200 CPO's, we need only compute 7 
of them. Let n denote the number of subjects in each group. Then the 
computation of the CPO for 5 = 1,...,50 is based on n (49,50,50,50) 
and X = (0, 2, 10, 30); the CPO for 5 = 51 , 52 is based on n = (50, 49, 50, 50) 
and X = (0,1,10,30), etc. 

The values of the CPO's for each prior choice with different models are 
given in Tables 4-6. Usually we prefer to select a diffuse prior for model 
selection and let the data do the talking. In our examples, the analysis is 
relatively insensitive to the prior choice. They all suggest that the model 
with k = 2 improves upon the model with k = 1. Although the pseudo- 
marginal likelihood factor for the model with A: = 3 is slightly larger than 
that of the smaller model, the improvement is very marginal. This suggests 
either model, A: = 2 or /: = 3, may be acceptable for the data set. In fact. 
Table 7, which provides the Gibbs approximation to the pseudo-marginal 
likelihood for k = 1,2,..., 9, shows that the pseudo-marginal likelihood 
is essentially fiat for 2 < A: < 9. Table 8 indicates that the value of the 
coefficient A, important for low-dose extrapolation, is not sensitive to the 
model for > 3. 
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Table 1: Gibbs Approximation to the Bayes Estimates when k = 1 





ocj = 2, 


7j = 20 


ttj = 1, 


7i = l 


ttj = 1, 


7j = .01 


A) (SD) 


.022 


(.011) 


.011 


(.008) 


.012 


(.008) 


95% C.I. 


(.010, 


.050) 


(.005, 


.034) 


(.006, 


.034) 


(SD) 


.161 


(.009) 


.175 


(.007) 


.175 


(.007) 


95% C.I. 


(.138, 


.174) 


(.157, 


.185) 


(.157, 


.186) 


i(.5) (SD) 


.077 


(.015) 


.084 


(.015) 


.084 


(.016) 


95% C.I. 


(.066, 


.083) 


(.075, 


.088) 


(.075, 


.089) 


i(l) (SD) 


.148 


(.024) 


.160 


(.025) 


.160 


(.025) 


95% C.I. 


(.129, 


.159) 


(.145, 


.169) 


(.145, 


.169) 


i(2) (SD) 


.274 


(.039) 


.294 


(.040) 


.294 


/ r\Ar\\ 

(.040) 


95% C.I. 


(.240, 


.293) 








.OUt7 ) 


i(3) (SD) 


.380 


(.049) 


.406 


(.050) 


.406 


(.050) 


95% C.I. 


(.337, 


.405) 


(.373, 


.425) 


(.373, 


.426) 


P(0), P(l) 


.021 


.167 


.011 


.170 


.012 


.170 


P(2), P(3) 


.290 


.395 


.302 


.414 


.304, 


.415 


P5i)(.025), RSD{.05) 


.022 


.184 


.079 


.228 


.076 


.225 


RSD{.075),RSD{.10) 


.350 


.521 


.381 


.538 


.377 


.534 


RSD{.125),RSD{.15) 


.697 


.877 


.699 


.865 


.695 


.861 
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Table 2: Gibbs Approximation to the Bayes Estimates when A; = 2 





aj =2,7j = 20 




1 ^\■^ 

aj = l,7j- = .01 


A (SD) 
95% C.I. 


.019 (.009) 
(.010, .041) 


.011 (.008) 
(.006, .033) 


.011 (.008) 
(.006, .034) 


fix (SD) 
95% C.I. 


.037 (.021) 
(.008, .088) 


.029 (.026) 
(.004, .100) 


.032 (.026) 
(.004, .100) 


k (SD) 
95% C.I. 


.063 (.011) 
(.039, .079) 


.068 (.012) 
(.035, .083) 


.067 (.012) 
(.035, .083) 


i(.5) (SD) 
95% C.I. 


.034 (.014) 
(.023, .053) 


.031 (.015) 
(.021, .058) 


.032 (.015) 
(.021, .057) 


i(l) (SD) 

95/0 CI. 


.096 (.020) 


.093 (.021) 


.094 (.022) 


i(2) (SD) 
95% C.I. 


.279 (.038) 
(.251, .296) 


.281 (.038) 
(.258, .299) 


.282 (.038) 
(.257, .300) 


i(3) (SD) 
95% C.I. 


.491 (.057) 
(.440, .527) 


.502 (.057) 
(.449, .535) 


.501 (.057) 
(.449, .532) 


P(0), P(l) 
P(2), P(3) 


.019 .113 
.293 .504 


.011 .103 
.290 .510 


.011 .105 
.291 .510 


RSD{mb),RSD{.Oh) 
RSD{.075\RSD{.10) 
RSD (.125), RSD {.lb) 


.139 .479 
.715 .910 
1.083 1.240 


.288 .582 
.798 .980 
1.142 1.290 


.275 .568 
.785 .968 
1.131 1.280 
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Table 3: Gibbs Approximation to the Bayes Estimates when /: = 3 





aj = 2,7j = 20 


ttj = l,7j = 1 


aj = l,7i = .01 


k (SD) 
95% C.I. 


.018 (.008) 
(.010, .041) 


.012 (.008) 
(.006, .034) 


.011 (.008) 
(.006, .033) 


/3i (SD) 
95% C.I. 


.030 (.017) 
(.008, .071) 


.025 (.019) 
(.004, .074) 


.025 (.020) 
(.004, .073) 


k (SD) 
95% C.I. 


.025 (.014) 
(.005, .057) 


.023 (.018) 
(.002, .066) 


.023 (.018) 
(.002, .067) 


k (SD) 
95% C.I. 


-.019 (.006) 
(.007, .029) 


.020 (.007) 
(.003, .031) 


.020 (.008) 
(.003, .031) 


i(.5) (SD) 
95% C.I. 


.023 (.014) 
(.012, .041) 


.021 (.014) 
(.008, .043) 


.021 (.014) 
(.007, .041) 


i(l) (SD) 
95% C.I. 


.071 (.020) 
(.049, .100) 


.066 (.022) 
(.041, .101) 


.066 (.022) 
(.040, .100) 


i(2) (SD) 
95% C.I. 


.268 (.037) 
(.239, .291) 


.260 (.038) 
(.231, .288) 


.260 (.038) 
(.232, .289) 


i(3) (SD) 
yo /o ^^.1. 


.560 (.059) 


.555 (.061) 
( 500 601^ 


.558 (.061) 
( 502 .601) 


P(0), P(l) 
P(2), P(3) 


.018 .088 
.282 .572 


.012 .077 
.269 .565 


.011 .076 
.270 .567 


KTiP(.025),KTi:>(.05) 

RSD{mh),RSD[.ld) 
R£D{.nh),RSD[.\h) 


.195 .624 
.888 1.091 
1.258 1.405 


.370 .741 
.983 1.173 
1.332 1.472 


.381 .748 
.989 1.177 
1.335 1.474 
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Table 4: Gibbs Approximation to the CPO for the Prior with Qj = 2 and 



7j = 20 for all j 





s 


1-50 51-52 53-100 101-110 111-150 151-180 181-200 


ln(nc.) 




c. 


.980 .163 .833 .285 .708 .389 .601 
(.017) ^.024) (.024) (.037) (.037) (.047) (.047) 


-78.242 


k = 2 




.982 .110 .887 .288 .706 .495 .494 
(.015) (.022) (.022) (.037) (.037) (.055) (.055) 


-72.743 


k = Z 




.982 .084 .912 .277 .717 .563 .424 
(.016) (.021) (.022) (.036) (.036) (.059) (.058) 


-70.863 



Table 5: Gibbs Approximation to the CPO for the Prior with = 1 and 



7j = 1 for all j 





s 


1-50 51-52 53-100 101-110 111-150 151-180 181-200 


ln(nc.) 


k = l 




.988 .166 .830 .297 .696 .406 .584 
(.014) (.025) (.025) (.039) (.039) (.049) (.049) 


-77.560 


k = 2 


SD{c.) 


.989 .100 .895 .285 .709 .501 .488 
(.015) (.022) (.023) (.038) (.038) (.056) (.056) 


-71.942 


k = ^ 


c. 


.988 .071 .923 .264 .730 .555 .431 
(.015) (.022) (.033) (.037) (.037) (.060) (.060) 


-70.132 
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Table 6: Gibbs Approximation to the CPO for the Prior with aj = 1 and 
7j = .01 for all j 





8 


1-50 51-52 53-100 101-110 111-150 151-180 181-200 


ln(nc.) 


k = l 


c, 

SD(c.) 


.988 .167 .829 .298 .695 .407 .582 
(.015) (.025) (.025) (.039) (.039) (.049) (.049) 


-77.623 


k = 2 


SD{c,) 


.989 .100 .894 .286 .708 .501 .487 
(.015) (.022) (.023) (.038) (.038) (.056) (.056) 


-71.997 


k = Z 


SD{c,) 


.988 .073 .922 .264 .730 .556 .431 
(.015) (.023) (.024) (.037) (.037) (.061) (.060) 


-70.122 


Table 7: Gibbs Approximation to the Pseudo-Marginal Likelihood In C{k) 


k 


123456789 


Qj = 2 & 7j = 20 

Qj = 7j = 1 
Q, = 1 t 7, = .01 


•78.242 -72.743 -70.863 -70.358 -70.368 -70.784 -71.334 -71.959 -72.912 
-77.500 -71.942 -70.132 -69.856 -69.441 -69.764 -70.065 -70.058 -70.423 
-77.623 -71.997 -70.122 -69.798 -69.743 -69.603 -69.864 -70.044 -70.425 j 


Table 8: Gibbs Approximation to the Bayes Estimates of 0i 




k 


123456789 




ttj = 2 & 7j = 20 

Qj = 7; = 1 
ttj = 1 & 7j = .01 


.161 .037 .030 .027 .026 .026 .025 .024 .024 
.175 .029 .025 .024 .024 .024 .023 .024 .025 
.175 .032 .025 .025 .024 .024 .024 .024 .025 
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6. Concluding Comments 

We have illustrated the usefulness of the Gibbs sampling method in the 
linearized multi-stage model. The results can also be extended to the hier- 
archical Bayesian setup where the hyperparameters aj and yj are random 
with appropriate distributions. The implementation of the Gibbs method 
can be carried out by methods similar to those in Gelfand £ind Smith (1990). 
When there is no conjugacy structure involved, the Metropolis algorithm 
(Metropolis, et al. 1953) can be used for sampling the variates. The re- 
sults can also be extended to different prior distributions; for example, 
instead of gamma distribution, we can assume that each /3j has a truncated 
normal distribution, because /3j > 0. Then we cein specify that the mean 
and variance of the normal distribution have normal and inverse gamma 
distributions. The Gibbs method can be implemented as in Section 2. 
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